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Abstract 

We perform molecular dynamics simulations of 512 water-like molecules that interact via the 
TIP5P potential and are confined between two smooth hydrophobic plates that are separated by 
1.10 nm. We find that the anomalous thermodynamic properties of water are shifted to lower 
temperatures relative to the bulk by ~ 40 K. The dynamics and structure of the confined water 
resemble bulk water at higher temperatures, consistent with the shift of thermodynamic anomalies 
to lower temperature. Due to this T shift, our confined water simulations (down to T = 220 K) 
do not reach sufficiently low temperature to observe a liquid-liquid phase transition found for bulk 
water at T ~ 215 K using the TIP5P potential. We find that the different crystalline structures 
that can form for two different separations of the plates, 0.7 nm and 1.10 nm, have no counterparts 
in the bulk system, and discuss the relevance to experiments on confined water. 
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I. INTRODUCTION 



the topic 



Despite the numerous accomplishments in water research to date 
continues to be the subject of intense interest. In particular, water confined in nanoscale 
geometries has garnered much recent attention due to its biological and technological impor- 
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tance |5|, I6JI • Confinement can lead to changes in both structural and dynamical properties 
caused by the interaction with a surface and/or a truncation of the bulk correlation length. 
Moreover these changes depend on whether the interactions of water with the wall particles 
are hydrophilic or hydrophobic 7]. One of the motivations for studying two different kinds of 
interactions arises from studies of protein folding, since the folding of a protein is influenced 
by its hydrophobic and hydrophilic interactions with water jjj. 

It is not clear exactly how the dynamics of liquids depends on the nature of the confining 
surfaces. The behavior may change depending on the surface morphology. Simulations 
of simple liquids show that the dynamics typically slow down near a non-attractive rough 
surface while the dynamics speed up near a non-attractive smooth surface M. A slowin 



down of water dynamics near a 



hydrophilic surface has been experimentally observed 



vmg 

s 



Water confined in Vycor ll|, [l2j has at least two different dynamical regimes arising from 



the slow dynamics of water near the surface and fast dynamics of water far away from the 



surfaces 



13, 



14. 
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One anomaly hypothesized to occur in supercooled water is the emergence of a phase 
transition line separating liquid states of different densities. This phenomenon is called a 
liquid-liquid (LL) phase transition 19, 2(1 EjE, 3- 24]. A LL transition has been seen in 



a variety of simulation models of water jl9j, |2JJ |26J] , but is difficult to observe experimentally 
due to the propensity of ice to nucleate at temperatures where a transition is expected. 
Nonetheless, indirect evidence of a transition has been found [2( 
some simple models of liquid also show a LL phase transition_ 

Bulk water simulations using the TIP5P potential 0, |3j| indicate the presence of a 
jL pha se transition ending in a second critical point at T rs 217 K and p ~ 1.17 g/cm 3 





25l l38| . A LL phase transition has been suggested based on simulations using the ST2 



potential confined between smooth plates 



simulations using the TIP4P potential 



theoretical work 



40 



39 



41 



. A liquid-to-amorphous transition is seen in 
42 1 confined in carbon nanotubes h!3j . Recent 



45j suggests that hydrophobic confinement suppresses the LL transition 
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to lower T. Here we aim to determine how confinement between smooth hydrophobic walls 
affects the location of the the LL critical point as well as the overall thermodynamic, dynamic 
and structural properties. 

The freezing of water in confined spaces is also interesting. On one hand, recent experi- 
mental studies of water confined in carbon nanopores show that water does not crystallize 
even when the temperature is cooled down to 77 K Q]. On the other hand, computer 
simulation studies show that models of water can crystallize into different crystalline forms 



when confined between surfaces 



44, 



47 



4£ . 14911 . For example, monolayer ice was found in 
simulations using the T1P5P mode] of water 4^ . By applying an electric field along lateral 
directions (directions perpendicular to the confinement direction) another crystalline struc 
ture for three molecular layers of water confined between two silica plates was found 



Also, bilayer hexagonal ice was found in simulations using the TIP4P model 46]. In general 
these simulations predict a variety of polymorphs in confined spaces, but the crystalline 
structures found have yet to be observed in experiments. 

This paper is organized as follows: In Sec. II, we provide details of our simulations and 
analysis methods. Simulation results for the liquid state are provided in Sees. Ill, IV, and 
V. The crystal states are discussed in Sec. VI, and we conclude with a brief summary in 
Sec. VII. 



II. SIMULATION AND ANALYSIS METHODS 

We perform molecular dynamics (MD) simulations of a system composed of iV = 512 
water-like molecules confined between two smooth walls. The molecules interact via the 



TIP5P pair potential 



which, like the ST2 [5fJ potential, treats each water molecule as a 



tetrahedral, rigid, and non-polarizable unit consisting of five point sites. Two positive point 
charges of charge qn = 0.241e (where e is the fundamental unit of charge) are located on 
each hydrogen atom at a distance 0.09572 nm from the oxygen atom; together they form an 
HOH angle of 104.52°. Two negative point charges (q e = —qn) representing the lone pair 
of electrons (e~) are located at a distance 0.07 nm from the oxygen atom. These negative 
point charges are located in a plane perpendicular to the HOH plane and form an e~Oe~ 
angle of cos -1 (l/3) = 109.47°, the tetrahedral angle. To prevent overlap of molecules, a fifth 
interaction site is located on the oxygen atom, and is represented by a Lennard- Jones (LJ) 
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potential with parameters oqo — 0.312 nm and too — 0.6694 kJ/mol. 



The TI 

is present 



5P potential accurately reproduces many water anomalies when no confinement 



251 ]. For example, it accurately reproduces the density anoma 



P = 1 atm. Its structural properties compare well with experiments 



25 



vat T = 277 K and 



36, 
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5lj. TIP4P 



and TIP5P are known to crystallize |25j, |52j within accessible computer simulation time 
scales; TIP5P shows a "nose-shaped" curve of temperature versus crystallization time j^J, 
a feature found in experimental data on water solutions 53]. TIP5P simulations also show 
a van der Waals loop in the P — p plane at the lowest T accessible with current computation 
facilities j^j]. This loop indicates the presence of a first-order LL transition. Ref. [3] 
estimates that a LL-transition line ends in a LL critical point C located at Ta = 217 ±3 K, 
P& = 340 ± 20 MPa, and p c > = 1.13 ± 0.04 g/cm 3 . 

In our simulation, water molecules are confined between two infinite smooth planar walls, 
as shown schematically in Fig. ^ The walls are located at z = ±0.55 nm, corresponding to 
a wall- wall separation of 1.1 nm, which results in « 2 — 3 layers of water molecules. Periodic 
boundary conditions are used in the x and y directions, parallel to the walls. 

The interactions between water molecules and the smooth walls are designed to mimic 
solid paraffin Sand are given b y 



U (Az) = 4e 



ow 



/ ^ow N 9 _ / cow V 
\Az) \ Az ) 



(1) 



54, 



56]. 



Here Az is the distance from the oxygen atom of a water molecule to the wall, while eow = 
1.25 kJ/mole and crow — 0.25 nm are potential parameters. The same parameter values 
were used in previous confined water simulations 

We perform simulations for 56 state points, corresponding to seven temperatures T = 220, 
230, 240, 250, 260, 280 i and 300 K, and eight densities p = 0.80, 0.88, 0.95, 1.02, 1.10, 1.17, 
1.25, and 1.32 g/cm 3 |57J. The range of density values takes into account the fact that the 
water-wall interactions prevent water molecules from accessing a space near the walls. Our 
determination of p is discussed in detail in the next section. The raw "geometric" densities 
used are p = 0.60, 0.655, 0.709, 0.764, 0.818, 0.873, 0.927 and 0.981 g/cm 3 . 

For each state point, we perform two independent simulations to improve the statistics. 
We control the temperature using the Berendsen thermostat with a time constant of 5 ps 



|58| and use a simulation time step of J. 
interactions we use a cutoff of 0.9 nm 



s, just as in the bulk system |25(. For long-range 



We calculate the lateral pressure P\\ = (P xx + P yy )/2 using the virial expression for 
the x and ^-directions 59]. We obtain the pressure alon g th e transverse direction, P± by 
calculating the total force F waU perpendicular to the wall 39], 



± 



wall | 



L X Ly L X Ly 



(2) 



Here, F i wa n is the force produced by oxygen atom of water molecule i on the wall. Hydrogen 
atoms do not interact with the wall. In agreement with the simulations of Ref. 56j] using 
the TIP4P model for water and the water- wall interaction given by Eq. (Q), we find that the 
hydrogen atoms of the water molecules near a wall tend to face away from the wall, forming 
bonds with other molecules. 



III. PROPERTIES OF TIP5P CONFINED WATER 
A. Transverse Density Profile 

One of the problems when dealing with liquids in a confined geometry is how to define the 
density in a consistent way. Using a geometric definition p = Nm/L x L y L z (where m is the 
water molecule mass) underestimates the effective density since the repulsive interactions 
with the walls prevent molecules from coming too close to the walls. Hence we want to 
quantify the effective distance L' z perpendicular to the walls accessible to the water molecules 
(Fig. 1), and thus obtain a definition for p which is more readily comparable with the density 
of a bulk system. To estimate L' z , we calculate the density profile p(z) defined as the density 
of oxygen centers at z, shown in Fig. |2] for different temperatures and densities. In all cases 
studied, we observe that molecules cannot access the total available space between the walls, 
and that the accessible space along the transverse direction does not strongly depend on T 
and p. Hence we estimate 

T , T ^OW + fOO ncm (i\ 

L z = L z = 0.819 nm (3) 

independent of T and p; this leads to the effective density 

Nm 

P = TTJT- ( 4 ) 

Figure Efa) shows the effect on p(z) of changing p at T = 230 K. Since the typical 
oxygen-oxygen separation for nearest neighbor in bulk water is 0.28 nm, for the effective 



wall separation L' z = 0.819 nm one would expect that at most three water layers can be 
accommodated between the walls. At p = 1.02 g/cm 3 and p = 1.25 g/cm 3 , p(z) shows three 
clear maxima indicating the presence of a trilayer liquid. The two maxima next to the walls 
are the result of water- wall interaction. As density decreases below p = 1.02 g/cm 3 , the 
central maximum becomes nearly uniform and, at density p = 0.95 g/cm 3 , only the two 
maxima located next to the walls remain. This density corresponds to the least structured 
liquid. Upon further expansion, the structure of the liquid starts to increase since the 
bilayer splits into two sublayers for the lowest density p = 0.88g/cm 3 . As we will see in the 
next section, T = 230 K is below the temperature of maximum density (Tmd)- Hence the 
increase in the structure upon expansion corresponds to the anomalous decrease in entropy 
upon expansion found in bulk water below the T MD , 



Figure 0(b) shows the effect on p(z) of changing p at higher T, T = 300 K. The densities 
in Fig. Efb) are the same as those shown in Fig. |2fa). Similar to findings at T = 230 K, (i) 
the liquid at T = 300 K and p = 1.25 g/cm 3 is characterized by a trilayer structure and (ii) 
reducing the density transforms the trilayer liquid into a bilayer liquid at p = 0.88 g/cm 3 . 
However, comparison of Fig. Efa) and fefrho-z(b) shows that, at the lowest p studied, the 
layers of the bilayer structured liquid at T = 300 K do not split into two sub-layers, as is 
the case of T = 230 K since T = 300 lies above the Tmd- In fact, we find that at T = 230 K 
the sublayers present at p = 0.88 g/cm 3 merge into a single layer at p = 0.95 g/cm 3 , and 
the resulting p(z) resembles that shown in Fig. |2fb) at p = 0.88 g/cm 3 . 

The effect of changing temperature at low density is shown in Fig. Efc). As discussed 
above, at low density and low T we observe a bilayer liquid where the two layers split into two 
sublayers. Increasing T smoothes features of the density profile, namely (i) the splitting of 
the sublayers disappears, and (ii) the minimum of p(z) at z = nm becomes less pronounced 
since increasing T increases the entropy. While a bilayer — > trilayer liquid crossover is found 
upon isothermal compression, isochoric heating does not have such an effect. 

Figure |2fd) shows the effect of changing T at the high density, p = 1.25 g/cm 3 . At high 
p, the molecules are not able to displace perpendicular to the walls and the density profile 
is almost T independent. Indeed, the liquid has a trilayer structure at all T studied. 





6 



B. Low Temperature Phase Diagram 



Next we test how confinement affects the location in the P — T phase diagram of the LL 
phase transition line and the second critical point found in bulk water simulations using the 



TIP5P model 



38| . To determine the phase behavior, we evaluate P as a function of p. 



If there is a second LL phase critical point, then it should manifest itself in a van der Waals 
"loop" along isothermal paths at low T. 

Unlike bulk liquid systems, one must be careful to interpret separately the results in- 
volving the lateral pressure P\\ and those involving the transverse pressure P±, since the 
thermodynamic averages of these quantities will usually be different. Phase separation will 
only be apparent in P, since the separation of the plates is too small to allow for the exis- 
tence of two distinct phases in the transverse direction. In Figs. Eta) and^b), we show P± 
and the lateral pressure Pi as functions of density along all seven isotherms. 

Since there can be no phase separation in the transverse direction, P± is a monotonically 
increasing function of the density (Fig. 0(a)). P\\ is also a monotonic increasing function of p, 
but as T decreases, the isotherms of Pi become "flatter" in the region near p w 0.95 g/cm 3 
(Fig. [2Kb)). The presence of this region is consistent with the possible existence of a van 
der Waals loop at lower T and, therefore, is consistent with the possible existence of a LL 
phase transition line ending at a second critical point at a value of T lower than 220 K, 
the lowest simulated temperature. At 220 K no phase separation occurs, consistent with 
the simulations of bulk water where the same TIP5P potential gives a second critical point 
with Tc = 217 ± 3 K [25]. While we are unable to simulate the temperatures below Tc 
of the bulk system, comparison of the lowest T isotherm with that of Ref. j3, El suggests 
that we will need to go well below the bulk Tc to see phase separation. Thus our results 
suggest that the presence of hydrophobic walls shifts a possible second critical point to lower 
T. Along an isothermal path, the critical point can be located by the point where the slope 
(dP/dp)x and curvature (d 2 P / dp 2 )x simultaneously equal zero. We can estimate this point 
by plotting the values of the minimum slopes along each isotherm and extrapolating the 
slopes to find the T at which the slope is zero. This estimate yields critical temperature 
T c > « 162 ±20 K (Fig. HI). 

Figures Etc) andOtd) show the T-dependence of P± and P for different densities. P±(T) 
is a monotonic function of T for all p. Similar behavior is observed for P\\(T) at large p. 
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However, for 0.88 g/cm 3 < p < 1.17 g/cm 3 , the isochores in the P\\ — T plane display minima, 
indicating the presence of a Tmd line, defined as the locus of points where (dP/dT)v = 
55| . For T < T MD , water confined between hydrophobic walls is anomalous, i.e., it becomes 
ess dense upon cooling. A Tmd line has also been found in TIP5P bulk water simulations 
25I . Comparison of Fig. EJc) and Fig. 2(a) of Ref. [2^ shows that the T MD locus in confined 



f. Q; 



water shifts to lower T. We also plot the Tmd for bulk water (from Ref. |25|) and confined 
water in Fig. A +40 K temperature shift in the Tmd of confined water overlaps these loci. 
Thus the effect of the hydrophobic walls in our system seems to be to shift the P — T phase 
diagram by AT —40 ±5 K with respect to bulk water. This is consistent with the second 
critical point shifting to lower T. 

Figure El shows the calculated potential energy for the lowest simulated temperature 
T = 220 K. We note two minima for the densities around p = 0.88g/cm 3 and p = 1.39g/cm 3 
respectively. Since the free energy F is given by 

F = K + U -TS, (6) 

where K, U, and S are the kinetic energy, potential energy and entropy respectively, at 
small T an extremum in U suggests an extremum in F. Hence the emergence of two minima 
at small T further supports the possibility of two stable liquids at low and high densities 
respectively. 

IV. STATIC STRUCTURE 

A. Radial Distribution Function 

In Sec. II we studied the structure of water along the direction perpendicular to the walls. 
To aid in comparing the structural properties with those of bulk water, we next focus on 
the lateral oxygen-oxygen radial distribution function (RDF) defined by 



d f k - Zj\ + Y)~ d y Zi ~ Zj][ ~ y 



(7) 



Here V is the volume, is the distance parallel to the walls between molecules % and j, Zi 
is the ^-coordinate of the oxygen atom of molecule i, and 5(x) is the Dirac delta function. 
The Heaviside functions, 9(x), restrict the sum to a pair of oxygen atoms of molecules 
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located in the same slab of thickness 8z = 0.1 nm. The physical interpretation of g\\(r) is 
that g\\(r)2iirdr5z is proportional to the probability of finding an oxygen atom in a slab of 
thickness 5z at a distance r parallel to the walls from a randomly chosen oxygen atom. In 
a bulk liquid, this would be identical to g(r), the standard RDF. 

Figure [7(a) shows the effect on g\\(r) of increasing p at low T . At low p, the bilayer liquid 
is characterized by a RDF that resembles the RDF of bulk water at p = 0.99 g/cm 3 [60], 
with maxima at r « 0.28, 0.45, and 0.67 nm. The well-defined maxima and minima in both 
p(z) (Fig. 12(a)) and g\\(r) (Fig. [7(a)) indicate that at low-T and low-p the liquid is highly 
structured, with a structure parallel to the walls similar to corresponding bulk liquid water 
at low density |6l|. As p increases, the liquid becomes less structured as indicated by the 
decreasing height of the second and third peaks of g\\(r), and by the disappearance of the 
sublayers in p(z) (Fig.[2(a)). A comparison of g\\(r) at p = 1.25 g/cm 3 with the g(r) for bulk 
water at high density from Ref. 6l( shows that the distributions are very different. However, 
we find that if we calculate g\\(r) for only —0.17 nm < z < 0.17 nm (which corresponds to 
the location of the central layer), then g\\(r) at low-p and high-p resemble g(r) for bulk water 
at both low and high densities respectively. Thus the evolution of the structure parallel to 
the walls of the central layer mimics the structural changes when going from low-density 
bulk water to high-density bulk water. 

The effect of increasing the density at T = 300 K is shown in Fig. [7(b). The main effect of 

compression is to "redistribute" the molecules parallel to the walls. At p = 0.88 g/cm 3 , g\\(r) 

is similar to the distribution shown in Fig. E(a) at the same density, which resembles g(r) 

for low-density bulk water. However, the maxima and minima in g» (r) are less pronounced 

at T = 300 K. At intermediate p, g\\(r) becomes very close to 1 for r > 0.45 nm. As density 

increases up to p = 1.25 g/cm 3 , a weak peak appears at r ~ 0.9 nm and a maximum occurs 

at 0.55 nm. Furthermore, the first maximum of g\\{r) becomes wider and the first minimum 

shift towards r « 0.4 nm. The resulting g\\ (r) at high p and high T has many features of g(r) 

| I 

for bulk water obtained experimentally at T = 268 K |6J|. For example, the oxygen-oxygen 
g{r) in Ref. j^l] shows a weak peak at r rs 0.9 nm, a clear maximum at 0.61 nm, and the 
first minimum is located at r « 0.41 nm. Furthermore, a shoulder in the first maximum of 
g(r) develops at high density which is consistent with the increase of the width of the first 
peak of g\\{r) in Fig. [7(b) as p increases. 

To complete the comparison of the structure of confined water with the structure of bulk 
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water, we also evaluate the T-dependence of g\\(r). In Fig. Efc), we show g\\{r) for various 
T at p = 0.88 g/cm 3 . The effects of increasing T at low p are similar to those observed in 
Fig. [7(a) when increasing p at low T. More specifically, the minima and maxima of g\\(r) 
become less pronounced as T increases and g\\(r) becomes flatter for r > 0.5 nm. We note 
that similar changes are also found in p(z) when (i) increasing T at low p (Fig. 12(a)) and 
(ii) increasing p at low T (Fig. 12(c)). 

Figure [7(d) shows the effect of T on g\\(r) when the density is fixed at p = 1.25 g/cm 3 . 
The effects of increasing T at high p (Fig. [7(d)) are similar to the effects of increasing p at 
a high T (Fig. E(b)). More specifically, the first minimum of <7||(r) shifts to a larger r and 
becomes less pronounced, while the second and third peaks located at 0.45 nm and 0.7 nm, 
respectively, merge and form an intermediate peak at r = 0.55 nm. This suggests that as p 
increases at high T, the preferred distance between second neighbors increases (parallel to 
the walls) and local tetrahedral order decreases. The emergence of a peak at 0.55 nm on 
heating at 1.25 g/cm 3 is in contrast to the behavior of bulk water, where the disappearance 
of the peaks at 0.45 and 0.7 nm gives rise to nearly featureless behavior of g(r) beyond the 
first peak (see, e.g., Ref. j3]). Hence at high T, confinement gives rise to structure that 
would not be present in bulk systems, presumably because molecules orient relative to the 
walls. It is interesting to note that, as shown in Sec. [21 there is almost no change in p(z) with 
T, indicating that the rearrangement of the molecules parallel to the walls has no effect, on 
average, on the organization of molecules perpendicular to the walls. 

B. Static Structure Factor 

Next we calculate the lateral static structure factor S\\ (q), defined as the Fourier transform 
of the lateral radial distribution function g\\ (r) , 

^teH^E^-^}. (8) 

*>i 

where the g-vector is the inverse space vector in the xy plane and r is the projection of 
the position vector on the xy plane. The structure factor will be particularly useful for 
comparison with the crystal structure in Sec. IV1[ where distinct Bragg peaks in S\\ (q) appear. 

Figures |H(a) and [S(b) show the effect of density on the lateral structure factor S\\ (q) 
for T = 230 K and T = 300 K. The structure factor of low-density and low-temperature 
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confined water is similar to bulk water. The presence of a "pre-peak" in Su (q) at low p can 
be attributed to the existence of pronounced tetrahedral order in the low-temperature liquid. 
At high p or high T, this feature is reduced, as the molecular order becomes less tetrahedral 
and core repulsion dominates. This behavior is similar to bulk water, but with a shift to 
lower densities. We show the evolution of Sn (q) as a function of T for two different densities 
p = 0.88 g/cm 3 and p = 1.25 g/cm 3 in Figs.|Sfc) and|Hfd). The structure of low-temperature 
water for p = 0.88 g/cm 3 is similar to low-density water. When the temperature is increased, 
the repulsive region of the potential begins to dominate and tetrahedrality is reduced. The 
first two peaks in the structure factor merge to form a single peak (Fig. Efc)). However, at 
high density p = 1.25 g/cm 3 , a change is temperature does not change the structure factor 
significantly (Fig. Et^d)). Similar behavior is seen in bulk water jf^]. 



V. DYNAMICS 



Thus far, we have seen that if a LL transition exists for confined water, it is shifted to 
lower T than for bulk water, and that the tetrahedral order that gives rise to this behavior is 
also suppressed. Hence it is natural to consider whether the dynamic properties of confined 
water exhibit the same temperature shift found for the thermodynamic properties relative 
to bulk water. For example, how is the maximum in diffusivity under pressure shifted under 
confinement? To compare with the bulk system, we calculate the lateral mean square dis- 
placement (MSD). We can evaluate the diffusion coefficient D from the asymptotic behavior 
of the lateral MSD using the Einstein relation 

(rf ) = 2dDt, (9) 

where (r|) is the mean square displacement parallel to the walls over a given time interval 
t, and d is the system dimension Since we calculate the diffusion only in the lateral 
directions, d = 2. 

Figures Efa) andEfb) show the dependence of the lateral MSD on p at fixed T = 230/^ 
and T = 300-fT. We also plot the dependence of the lateral MSD on T for two different tem- 
peratures in Figs. Etc) and^d), using a log-log scale to emphasize the different mechanisms 
seen on different time scales: 

(i) An initial ballistic motion, where the lateral MSD is a quadratic function of time, 
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W(*)> ~ t\ 



(ii) An intermediate "flattening" of the lateral MSD, due to the transient caging of 
molecules by their hydrogen bonded neighbors. This effect is most noticeable at the 
lowest T studied, and does not occur at high T. 

(iii) Long time scales on which particles diffuse randomly, and so (r|(i)) ~ t. 

To determine whether there is an anomaly in the density dependence of D, we plot D 
along isotherms, in Fig. HUT a) for confined water and in Fig. ITUT b) for bulk water. For 
T < 250 K, we find that D has a maximum at p pa 1.05 g/cm 3 . In bulk water, a similar 
behavior is found, but at T pa 290 K, 40 K higher than the confined system. Moreover, 
this shift in a dynamic anomaly is consistent with the shift of thermodynamic anomalies. 
Qualitatively, the maximum in D can be understood as a competition between weakening or 
breaking of hydrogen bonds under pressure jf^ (which increases D) and increased packing 
(which reduces D). 

We next study the effect of confinement on the rotational motion of water molecules. The 
rotational motion was analyzed by calculating the rotational autocorrelation time for all the 
state points and are compared with the rotational autocorrelation time for bulk water for 
few state points. The rotational autocorrelation function C±(t) is defined as 



where e*j(t) is the unit dipole vector of molecule % at time t. For large times, C\(t) can be 
fit with a stretched exponential function 



where A, [3 < 1 are constants and r = r(p, T) is the orientational autocorrelation time, 
which depends on both density and temperature. In Fig. ^2 we show Ci(i) for different 
temperatures and densities. Figure 12(a) shows the inverse of the orientational autocorrela- 
tion time which is proportional to the rotational diffusion. For a comparison with the bulk 
water rotational diffusion, r _1 is also shown in Fig. 12(b). Both the translational diffusion 
constant [Fig. ITUT a)] and the inverse of the orientational autocorrelation time [Fig. 12(a)] 
show similar behavior. Similar results have been found for bulk water with the SPC/E 




(10) 



C x {t) = Ae^ 13 



(11) 
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potential For low temperature, the maxima occur at the same density for D and r . 
However at high temperatures (T > 260K), where the D is a monotonically decreasing 
function of density, r -1 has a maximum and a minimum similar to its bulk counterpart 
[Fig. 12(b)]. 

VI. CRYSTALLIZATION OF TIP5P CONFINED WATER 

Bulk TIP5P water crystallizes within the simulation time at higher densities and, for a 
given density, the crystallization time has a minimum at T = 240 K 25j] . We next investigate 
whether crystallization occurs in confinement, and whether the structure differs due to 
the surface effects. It has been found experimentally that water confined in hydrophobic 
carbon nanopores does not crystallize, even at very low temperatures lo| . However, the 
crystallization of confined water is seen in some simulations 0, 47]. We find that our 
system crystallizes to what appears to be a trilayered ice structure at high density and 
that the resultant ice has a density p = 1.32 g/cm 3 . A similar crystallization appears in 
simulations when an electric field is applied in a lateral direction |48j|. The surface in this 
kind of simulation exhibits the embedded crystal structure of silica, where the oxygen and 
silicon atoms are arranged in out-of-registry order. This suggests that the crystalline form 
we find in confinement does not depend on the morphology of the surface. We show the 
structure of the ice and the static lateral structure factor in Fig. ED 

To investigate whether confined water crystallizes in the same crystalline form when the 
separation between the plates is different, we repeat our simulation for a plate separation 
of 0.7 nm, and again the system crystallizes at T = 260 K. The system crystallizes into a 
monolayer ice, also seen in simulations in Ref. j^|. The ice structure and its lateral static 
structure factor are shown in Fig. The density of monolayer ice is p = 0.93 g/cm 3 . We 
list the temperature, pressure, and potential energy for these crystals in Table I. 

VII. CONCLUSIONS 

We have systematically investigated the effect of confinement on TIP5P water between 
two parallel smooth hydrophobic plates, separated by 1.10 nm. We found that the overall 
phase-diagram is shifted to lower temperature and lower density compared to bulk TIP5P 
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water by AT rs — 40K. The shift to lower temperature compared to the bulk water can be 
understood qualitatively. Since the confinement walls do not form any hydrogen bonds with 
water molecules the average number of hydrogen bonds per molecule in confined system is 
smaller than in bulk water. This is analogous to having bulk water at high temperatures. 

We do not see a LL phase transition for the state points we have been able to simulate, 
but we do see a pronounced inflection in the P — p isotherms (Fig. |3J), which is consistent 
with a LL phase transition at lower T. Since the phase diagram is shifted ~ 40 K lower 
in temperature, our results are consistent with the possibility that there indeed is a LL 
transition at a temperature too low to simulate. This shift of thermodynamics qualitatively 



agrees with the theoretical predictions of Ref. 45]. The structure of confined water is similar 



to the structure bulk water at a lower density, and shows a similar evolution of the structure 
with changes in density and temperature. At a given temperature, as the density increases, 
water changes from a bilayer liquid (at low density) to a trilayer liquid (at high density). We 
find that the confinement affects the translational diffusion as well as the rotational motion 
of water molecules. The rotational diffusion anomaly precedes the translational diffusion 
anomaly, just as occurs for bulk water. 

We were able to crystallize water for a few state points. It crystallizes spontaneously to 
a trilayer ice at T=260K. Monolayer ice was formed when the separation between the plates 
was decreased to 0.7 nm. The crystalline structures are different from the polymorphs of 
bulk water, and should be relevant for confined water. 
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TABLE I: Thermodynamic properties of the crystals formed at T = 260K 
Plate Separation (nm) Pressure (MPa) Potential Energy (kJ/mol) 
0.7 -90.30 -39.86 

1.1 652.89 -44.08 
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Water MdecuhK 

FIG. 1: Schematic of the simulated system. Water molecules are confined between two smooth 
hydrophobic plates in an xy plane, separated by L z = 1.1 nm. The figure also indicates L' z , the 
effective length of confinment of the water molecules. 
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FIG. 2: Density of water p(z) along the confinement direction for different densities at constant 
temperatures (a) T = 230 K and (b) T = 300 K. (c) Density of water along the confinement direc- 
tion for different temperatures at constant densities (c) p = 0.88 g/cm 3 and (d) p = 1.25 g/cm 3 . 
See discussion in text. 
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FIG. 3: (a) Transverse pressure P± (perpendicular to the walls) and (b) lateral pressure Py (parallel 
to the walls) as functions of density p for all simulated T. For clarity, each curve is shifted by 
n x 50MPa for Pi for and n x lOOMPa for Py. While Pi shows no infection for all T, Py starts to 
flatten at p « 0.87 g/cm 3 for the lowest T presented. This is consistent with the possibility of a 
second critical point at lower than T = 220 K. (c) Pi and (d) Py as a function of temperature for 
different densities. Only Pi displays minima, indicating the presence of a Tmd- 
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FIG. 4: Minimum values of (dP\\/dp)T as a function of temperature obtained from Fig. 0(d). 
The critical temperature Tc of the critical point of the LL transition is crudely estimated by 
extrapolation. 
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FIG. 5: A comparison of Tmd line for bulk and confined water, (a) Isochores corresponding to the 
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FIG. 6: Potential energy as a function of p for the lowest temperature we simulate, indicating the 
possibility of two stable liquid states at p = 0.88g/cm 3 and p = 1.39g/cm 3 . The two minima are 
indicated by vertical arrows. 
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FIG. 7: Lateral oxygen-oxygen radial distribution function <7|| (r) for different densities at a constant 
temperature (a) T = 230K and (b) T = 300 K; g»{r) for different temperatures at a constant density 
(c) p = 0.88 g/cm 3 and (d) p = 1.25 g/cm 3 . 



25 



2.5 



1.5 



CT 

c/T 



0.5 



1 1 1 , 1 , , 1 1 1 1 

; (a) 

'_ f- 

n T = 230K 

i ■ /I 
• i / 1 


, , , | i i 




p=u.oog/cm 
-■ p=1.02g/cm 3 
- p=1.25g/cm 3 




■ \ V l 




- 








?"=r , , . i . . , , i 


. . . i , , 






T=230K 
T=250K 
T=300K 



q(nm ) 




100 

q (nm"') 



200 
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FIG. 9: Lateral mean square displacement (r^(t)) for different densities at constant temperatures 
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FIG. 10: Diffusion constant D as a function of density for various temperatures in (a) confined 
water and (b) bulk water. Confined water shows diffusion anomalies at lower temperatures than 
in bulk water. The first diffusion anomaly (the increase of diffusion with increase in density) in 
the confined system begins to appear when the temperature drops to T « 250 K, which is « 40 K 
lower than in bulk water (which displays diffusion anomalies at T s=s 290K). 
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FIG. 11: Dipole-dipole autocorrelation function C\{t) for different densities at constant tempera- 
ture (a) T = 230 K and (b) T = 300 K. Also shown is C\(f) for different temperatures at constant 
density (c) p = 0.88 g/cm 3 and (d) p = 1.25 g/cm 3 . 
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FIG. 12: Inverse orientational autocorrelation time as a function of density for (a) confined water 
and (b) bulk water. Inverse autocorrelation time for confined water follows a similar behavior as 
bulk water, but temperature shifted by « 40 K to lower T. 
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FIG. 14: Monolayer ice with L z = 0.7 nm and p = 0.93 g/cm 3 spontaneously formed at T = 260 K. 
(a) A lateral snapshot and (b) top view, (c) Lateral structure factor. 
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